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Abstract. Particle acceleration by means of non-linear plasma wave interactions 
is of great topical interest. Accordingly, in this paper we focus on the 
electron surfing process. Self-consistent kinetic simulations, using both relativistic 
Vlasov and PIC (Particle In Cell) approaches, show here that electrons can be 
accelerated to highly relativistic energies (up to lOOrrteC^) if the phase speed of the 
electrostatic wave is mildly relativistic (0.6c to 0.9c for the magnetic field strengths 
considered). The acceleration is strong because of relativistic stabilisation of 
the nonlinearly saturated electrostatic wave, seen in both relativistic Vlasov and 
PIC simulations. An inverse power law momentum distribution can arise for the 
most strongly accelerated electrons. These results are of relevance to observed 
rapid changes in the radio synchrotron emission intensities from microquasars, 
gamma ray bursts and other astrophysical objects that require rapid acceleration 
mechanisms for electrons. 
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1. Introduction 

Particle acceleration by means of nonlinear plasma wave interactions can arise 
from a wide variety of processes. These include wakefield accelerators; electron 
interactions with lower-hybrid wave packets and solitons; electron heating by 
collapsing electrostatic waves; and gyroresonant surfing acceleration P |2 E El ■ 
Here we focus on electron surfing acceleration (ESA). The principle of ESA is the 
trapping of electrons by a strong monochromatic electrostatic wave combined with the 
transport of the trapped electrons across an ambient magnetic field that is oriented 
orthogonally to the wavevector. As long as the electrons remain trapped in the 
potential of the wave, they are continously accelerated, orthogonally to the magnetic 
field direction and to the wave propagation direction The resulting phase space 

structures are thus not equilibrium structures, as in an unmagnetized plasma The 
electrons can be accelerated from the thermal population to energies at which they 
emit radio synchrotron radiation, in considerably less than an inverse ion gyroperiod 
|10|. The rapid electron acceleration makes ESA a candidate mechanism for particle 
acceleration at microquasars, for example GRS 1915+105 IT and gamma ray bursts 

m 

The electrostatic wave that traps the electrons can be excited by ion beams that 
arise ahead of shocks ^1 1141 1151 1161 [T7j . As the shock expands into the upstream 
plasma, it reflects a substantial fraction of the upstream ions. If the shock normal is 
oriented perpendicular to the upstream magnetic field, ions can be specularly reflected 
^lEl or they can undergo acceleration by shock surfing 118111911^ and other plasma 
processes |221 ■ The interplay between the plasma shock, the shock- refiected beam 
of ions and the electrons undergoing ESA has previously been analysed by means of 
kinetic particle-in-cell (PIC) simulations |23 ElESl ES EH ■ 

MHD shocks in the accretion disks of microquasars and in the outflow of gamma 
ray bursts may expand at a significant fraction of the speed of light c into the upstream 
region, as MHD simulations show for microquasars |28| and as experimental evidence 
indicates for gamma ray bursts The shock-reflected ions will move with, at least, 
that speed. Electrostatic waves driven by such ion beams have mildly relativistic 
phase speeds, which can increase their nonlinear stability and their lifetime compared 
to electrostatic waves with large but nonrelativistic phase speeds ;5i i29i i30i [3T) . 

In this paper, we examine the efficiency of ESA for waves driven by moderately 
dense proton beams that move at the mildly relativistic speeds Vb = 0.6c and Vh = 0.9c. 
The efficiency with which the electrons are accelerated by their cross-field transport 
will depend critically on the wave stability; the maximum electron energy will be 
between « 10 keV found for a proton beam speed of Vb ~ 0.06c IB^ and the GeV 
energies found for vi, — 0.99c JOj. We show below that electrons are accelerated up 
to Lorentz factors of the order of ten by ESA for Vb — 0.6c. The peak acceleration is 
further increased due to a strongly nonlinear density modulation of one proton beam 
that is associated with a jump in the beam velocity; this may evolve by a similar 
process to that reported in Ref. j33j for a Q-machine experiment and for the expansion 
of plasma into a vacuum j34' . This charge accumulation yields strong localized electric 
fields that can accelerate electrons up to 7 w 20, where the Lorentz factor is defined 
as 7(w) = (1 — v^/c^) with v being the electron velocity. The waves driven by 
the beam with Vb = 0.9c can accelerate electrons to 7 ~ 100, which is well above the 
Lorentz factors of the fast transient jets of the microquasar GRS 1915-1-105 and 
could also account for strong synchrotron emissions. In this regime of beam speeds, 
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the peak energy which the electrons can reach by ESA increases more than the beam 
speed. 

The calculation of wave stability, and thus the efhciency by which electrons can 
be accelerated, should be invariant with respect to the simulation scheme. However, 
the schemes that represent the plasma vary between different codes and this can 
sometimes lead to differences in the computed wave stability, e.g. due to different 
levels of simulation noise and perturbations of the islands of trapped particles 
The resulting different lifetimes may, in turn, affect the computed ESA efficiency. 
We, therefore, compare in section 2 the development of the beam instability in a 
PIC simulation with that in an electrostatic Vlasov simulation |31 OH for an 
unmagnetized plasma composed of ions and electrons. The PIC code advances a 
number of macro-particles along the characteristics of motion, given by 

at •jnie 

,2) 



for the electrons, and by 
dx p 



(3) 



dt ^nii 

for the protons. The Vlasov code solves the one dimensional relativistic Vlasov-Poisson 
system directly. This fully nonlinear self consistent system is governed by the Vlasov 
equation for the electron distribution function fi,{x,p,t) 

^ + J!—^_eE—=0 (5) 
dt me"f dx dp 

the Vlasov equation for the ion distribution function fi {x,p, t) 

'4 + ^f+eEf=0, (6) 
ot rui^ ox op 

and the Poisson equation for the electric field 

dE _ e 

dx Co 

Comparison of the two codes confirms almost identical nonlinear evolution of the 
electrostatic waves. We then perform electromagnetic PIC simulations with different 
strengths of the magnetic field in section 3, and discuss the results in section 4. 



.f,dv - / Jdv . (7) 



2. Simulation setup and code comparison 

As a perpendicular shock propagates into the upstream medium, beams are formed by 
the shock reflected protons. Ideally one would simulate the formation of these beams, 
but to approach this complete problem is, at the present time, too computationally 
demanding. Thus, we focus instead on a one-dimensional periodic system, in which 
the counter-propagating beams are present aXt — Q, the evolution of which is governed 
by the Buneman-type instabilities |13[I14) . Therefore, in this work we do not consider 
issues like the wave coherency orthogonal to the wavevector or effects introduced by 
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Figure 1. Schematic of the initial plasma distribution in our simulation. We 
model two beams that have an equal mean speed and move in opposite directions. 
The thermal spread of the background electrons and protons and that of both 
beams is small compared to the beam speed modulus. Beam 2 is hotter to account 
for the proton beam interaction with the upstream medium (the thermal spread 
in the simulations is less than shown here.). The initial mean speeds along Vy 
and Vz are zero for all particle species. The total proton density is equal to the 
electron density. 



the beam front, which arc the subject of future investigations. Referring to Fig. ^ 
beam 1 is composed of the protons that have just been reflected. This beam moves into 
the upstream plasma, is rotated by the upstream magnetic field B and returns to the 
shock as beam 2 [X7|. The background electrons and protons represent the upstream 
plasma prior to its shock encounter. In total, we thus have four plasma species, which 
are all modelled in a fully kinetic treatment due to the importance of the ion reaction 
to the wave fields 31 . We place our simulation box into the upstream plasma close to 
the shock. We neglect the fact that the reflection of protons due to the shock motion 
cannot produce perfectly anti-parallel beams, and align the mean velocity vectors 
of both beams with the simulation or s-direction. We also assume the plasma to be 
homogenous orthogonal to the beams. The mean speeds of the beams in the a;-dircction 
are vt for beam 1 and —vt for beam 2, so that the total current in the simulation box 
vanishes. Therefore, no return current develops for this initial condition. Initially, 
all components of the mean speed are zero for the background electrons and protons, 
as are the mean speed components of the proton beams orthogonal to the simulation 
box. We limit our simulation time to much less than an inverse proton gyrofrequency, 
thereby eliminating the effects arising from the magnetic rotation of the beam. The 
spatial scale of our simulation box is small compared to the forcshock region, so that 
we may take all four species as spatially homogeneous. 

As initial parameters for all simulations, we take the electron plasma frequency 

1/2 

i^p,e = ('T-ee^/eo™e) = 27r X 10^s~^. Here n^, e, eo and rUe denote the electron 
number density, the magnitude of the elementary charge, the dielectric constant and 
the electron mass, respectively. The number densities of beam 1 and beam 2 are 
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nti — ni,2 — rib = n^/W in the rest frame of the respective beams. Such densities 
for the shock reflected ion beams are representative for perpendicular shocks [67\ . 
The density of the background protons is Up = — 2j(vii)ni,, such that the net 
charge in the simulation box vanishes. The initial thermal spread of the electrons is 

For the simulations of beam instabilities in an unmagnetized plasma in this 
section, the background protons and beam 1 have the thermal speed vth.p = Vth,bi = 
VWvth,eirne/mp)'^^'^ , where rup is the proton mass. In all simulations, we use 
rrip/me — 1836. Proton beam 2 has the thermal speed Vth,b2 — 10wth,bi- The 
proton temperature is increased relative to the electron temperature to improve the 
resolution of the proton velocity distributions by the Vlasov code. For the simulations 
in a magnetized plasma, to be discussed in the next section, we use the values 

1/2 

Vth,p = vth,bi = VthA'^e/rrip) and Vth,b2 = lOwt/i.w. 

All particle species have a thermal speed that is small compared to the considered 
beam speeds vi, = 0.6c (slow beams) and Wf, = 0.9c (fast beams), and are thus cold. 
The temperatures of the background electrons and protons are higher than a few eV, 
and their number densities are lower than w 10^°to~'^, estimated for the accretion disk 
of an AGN . Our choice of initial conditions is due to restrictions in the computer 
time. No accurate phase space distribution functions are available for the beams at 
relativistic shocks and we thus take them to be cool, as in Ref. |^. In this section, 
our plasma is unmagnetized so that the developing instability is that discussed by 
Buneman in Ref. ^|, except for a factor of 7 (w;,) = 1/(1 — w^/c^), which has been 
examined numerically in Ref. |14| . To a good approximation, the most unstable 
Buneman wave is expected to grow at the wavenumber fc„ = ujp i^/v},, at a frequency 
u>u = i^p,e, and with the linear growth rate 

1 /3 

ujim ^ J^^{vb)i3V3ujl i,ujp,i/16) , (8) 

where Up^f, is the beam plasma frequency. 

We use M = 512 grid cells along the s-direction for both the Vlasov and PIC 
simulations. The total box length is L — 2A„ — Air/ku, and we normalize time by 
Tp = 27r/ujp^e- We analyse the electric field calculated by the simulations by first 
taking a Fourier transform over space and time. We separate the w, k spectrum into 
two sets of two quadrants that have either uj/k < Oorw/fc > 0. The consecutive 
inverse Fourier transforms applied to both data sets reconstruct the wave fields for 
the positive and negative phase speeds separately. In the absence of nonlinear mode 
coupling, we would thus separate the wave fields driven by both proton beams. We 
obtain the field amplitudes as a function of A:, i by applying 

M 

^^Ex{xi,t)exp{-ikjXi) (9) 



E{kj,t) = M 



-1 



to the electric field. The amplitudes of the electrostatic Buneman waves are given in 
physical units V/m. 



2.1. Slow beam 

The PIC simulation for Vb = 0.6c uses 1600 particles per cell (ppc) for the electrons and 
784 ppc for each of the proton species. The Vlasov simulation uses 8192 equidistant 
cells in px covering the momentum range —4.4 < (px/mec) < 4.4. For both 
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Figure 2. Comparison of waves driven by proton beams with vi, = 0.6c. The 
logarithmic amphtudcs of the waves with a negative phase speed {ku < 0) are 
shown in (a) and with fcu > in (b). PIC results are shown in blue and Vlasov 
results in red. 



simulations we calculate the amplitudes E(ku,t) and E{—ku,t) from the computed 
Ex using equation 1^1 and compare the results in figure |21 

The electric field in both simulations grows initially exponentially. The 
experimentally measured growth rate in the PIC simulation is Tex ~ 0.97LUim, with 
LUim/^p.e = 0.0228 inferred from Eq. |H1 In the Vlasov simulation the measured 
growth rate is practically identical. We observe from figure El that both waves in the 
Vlasov simulation grow symmetrically, whereas the wave with LUu/ku < in the PIC 
simulation grows earlier, and to a larger peak amplitude. 

Vlasov codes are free of noise, so we introduce an initial perturbation to encourage 
the growth of instabilities. We apply the initial perturbation to both beams in the 
Vlasov simulation at as in and therefore the initial electric field amplitudes for 
both waves are identical. The thermal spread of both beams and their difference is 
small compared to vt, and the growth rates of the instabilities driven by both beams 
are thus similar. The two waves in the Vlasov simulation thus grow at the same rate 
to a comparable saturation amplitude, as shown in figure [3 

In the PIC simulation, on the other hand, the initial wave amplitudes are given 
by noise, and the noise levels are higher for the hotter beam |22,. Therefore, the 
wave with uj/k < grows first. Since the weak wave with positive Uu/ku in the PIC 
simulation can interact nonlinearly with a limited phase space interval compared to 
that of the waves in the Vlasov simulation, the wave with negative oj^/k^ can grow 
to a larger amplitude without interacting with the second wave, e.g. by the formation 
of stochastic bands due to the simultaneous interaction of electrons with both waves 

SOI. 

The electric fields in both simulations saturate by the trapping of electrons as 
discussed for almost identical parameters in Ref. |29| and where the term trapping 
is used to denote electron reflection by the wave potential, as in As a result, 

phase space holes form in the electron distribution ^2 1421 I48| . The islands 
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Figure 3. Logarithmic plots of the electrostatic field energy density and the 
fc-spectrum for vi, = 0.6c. (a) The total electrostatic field energy density, (b) 
The power spectrum, integrated from ti to t2, of the total electrostatic field as a 
function of k/k^. Results from the PIC simulation are given in blue and from the 
Vlasov simulation in red. 

of trapped electrons have only a finite lifetime at this Vb |B3 ^^'^ collapse by the 
sideband instability [441 145| . Both simulations show a comparable lifetime for the 
nonlinearly saturated electrostatic wave. After t « 70Tp, the electric field amplitude 
fluctuates around a low and constant mean value. Further evidence for similar 
development of the electrostatic instability in both codes is provided by the energy 
density of the electrostatic field Etotai — M^^ X^jli ^^i^jji) the wave spectrum 
P{k) = {t2 - tiY^ E'^{k, t)dt, with ti = 200Tp and ^2 = 22QTp, which are shown 
in figure 13 The energy density of the electric field peaks at t « 30Tp and decreases 
thereafter, primarily because of the collapse of the electrostatic waves. At the end of 
the simulation, both codes show a comparable fc-spectrum at low wavenumbers. The 
PIC code gives stronger wave power at higher wave numbers, primarily due to its 
higher noise levels. 

The nonlinear interaction of the electrons with the electrostatic waves has 
increased the electron energy. Both simulation codes show, at t = 220Tp, a similar 
electron momentum distribution for px > 0, see Fig0] For Px < the electron 
momentum distribution calculated by the PIC code shows a tail. The higher peak 
momentum at negative Px in the PIC simulation is due to the larger amplitude of the 
wave with k = —ku compared to that in the Vlasov simulation, as shown in figure |21 
The larger amplitude leads to an island of trapped electrons that is wider in px] the 
collapse of the stronger wave thus redistributes the electrons over a wider px interval. 
The weaker wave with fc = /c„ in the PIC simulation cannot accelerate electrons to 
the same peak momentum as the waves in the Vlasov simulation, because initially it 
only accesses a smaller px interval. Conversely, the equally strong electrostatic waves 
with k = ku and fc = — fc„ in the Vlasov simulation lead to the observed symmetric 
momentum distribution. 
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Figure 4. The electron momentum distribution for vi, = 0.6c at t = 220Tp in the 
PIC simulation (blue) and in the Vlasov simulation (red). 

2.2. Fast beam 

The PIC simulation for vi, — 0.9c uses 3200 particles per cell (ppc) for the electrons and 
2048 ppc for each of the proton species. The Vlasov simulation uses 8192 equidistant 
cells in p^. covering the momentum range from —20 < (p^/m^c) < 20. We show the 
amplitudes E{ku,t) and E{ku/2,t) in figure|Sl 

The anticipated growth rate ujim ~ 0.0152, inferred from Eq. |H1 is accurately 
reproduced by the PIC simulation. In the Vlasov simulation, the maximum growth 
rate is Tex ~ 0.78wmi. From figure we find that initially the wave with k = ku, 
which is driven initially by the Buneman instability, is dominant, but at t ^ 200rp 
its amplitude decreases and a wave grows at ku/2. The coupling of the wave energy 
towards lower wavenumbers is a characteristic of the sideband instability |45j , and the 
decrease in the amplitude of the k — ku wave is associated with the collapse or merger 
of islands of trapped electrons, as discussed in Refs. |29[l4ti| . The time lag between the 
initial nonlinear saturation and the coupling of energy towards the sideband modes 
thus provides a measure of the non-linear stability of the electrostatic instability (which 
in the initial, linear phase is governed by the Buneman waves). Both codes show here 
a lifetime of about 150Tp. Figure |B1 shows that the energy density in the electrostatic 
field as a function of time is almost identical in both codes; so is the wave power, 
integrated between ti — 380Tp and t2 = 400Tp, at low wavenumbers. Again, the 
PIC code shows higher noise levels at larger wavenumbers. The electron momentum 
distributions at the end of the PIC and Vlasov simulations agree reasonably well, 
as shown in figure \7\ The electrons in both simulations reach a peak momentum 
\px\ ~ WnieC, however the Vlasov simulation has a flatter central maximum. The PIC 
code also evolves to a flat-top distribution, but more slowly. 
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Figure 5. Waves driven by the beams with vi, = 0.9c. The blue traces correspond 
to the waves with k = ku and the red to the waves with k = ku/'i. The ampHtudes 
of the waves with positive phase speed (fc > 0) are shown in (a) for the PIC 
simulation and in (b) for the Vlasov simulation. The amplitudes of the waves 
with negative phase speed (fc < 0) are shown in (c) for the PIC simulation and in 
(d) for the Vlasov simulation. 
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Figure 6. Logarithmic plots of the electrostatic field energy density and k- 
spectrum for ti;, = 0.9c in the PIC simulation (blue) and in the Vlasov simulation 
(red), (a) The total electrostatic field energy density (b) The power spectrum, 
integrated from 380Tp to 400Tp, of the total electrostatic field as a function of 
k jku- 



Efficient electron surfing acceleration 



10 



10° 




10 



i 



-20 



-10 



,0 

p / m c 



10 



20 



Figure 7. The electron momentum distribution for vi, = 0.9c at t = 400Tp in the 
PIC simulation (blue) and the Vlasov simulation (red). 



2.3. Particle distributions at late-time 

We note that the PIC work in Ref. yields a flat-top distribution at the later time 
t ~ lOOOTp, whose plateau has a width « GrTigC and which shows an electron phase 
space distribution similar to that of the Vlasov simulation shown in Fig. |21 The proton 
beam distributions computed by the Vlasov simulation, for the case of vt — 0.6c display 
ion holes which resemble qualitatively those in Ref. [221 • However, in the case of the 
Vlasov simulations presented here, ion holes develop almost symmetrically with both 
positive and negative phase speeds, rather than the 'one sided' distribution observed 
in PIC simulations. This is shown in Fig|Hl3 for the bulk protons and is attributable 
to the comparable amplitudes for instabilities associated with both beams. The two 
strong electrostatic waves can each decay into a low frequency mode with k — 2/c„ 
and a second Langmuir wave. In the PIC case, the ESWs associated with beam 2 
develop first, as discussed earlier. In both the slow and fast beam cases, and for 
both Vlasov and PIC codes, the proton distributions at late time are non thermal. In 
the absence of a magnetic field, their continuing evolution will be governed primarily 
by the effective collisionality, which naturally results from the finite resolution of the 
numerical schemes [21 BZ| , from particle collisions with the fluctuating electric fields 
0H1 and the stochastic interaction of the electrons with the electric fields of the ion 
beam structures |40| . Hence, the time asymptotic equilibrium shown in figure 8 is 
similar to the structural dissipative equilibria described in the Refs. pillIT) . 

3. Magnetic field effects 

In the preceeding section, we established that the beam-excited electrostatic waves 
can remain non-linearly stable for one hundred plasma oscillation periods or more. 
We have shown that the results from PIC and Vlasov simulations agree on this, which 
is evidence for a physical origin behind the wave stabilization. In the present section. 
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Figure 8. The final electron (a) and bulk proton (b) phase space distributions 
for Vf, = 0.6c at t 200Tp in the Vlasov simulation. The colour axis denotes the 
logarithmic magnitude of the distribution functions, /e and fi. 



we assess the efBciency of these waves with respect to trapping and ESA: this requires 
the inclusion of an external magnetic field. Under these conditions, our computational 
approach is based on PIC simulations. We use 512 simulation cells to resolve the x- 
direction of the simulation box with its length L = 47r/fc„. The minimum simulation 
wavenumber 2n/L < /c„ allows for the development of the sideband instability and 
the merger of trapped particle islands, and thus corresponds to the case of the long 
simulation box in Ref. I32j. The magnetic field orientation is aligned with the z- 
axis and this orientation, together with the relativistic wave speeds, makes the results 
directly comparable to the work in Ref. 0. 

We choose values for the magnetic field strength that correspond to an electron 
gyrofrequency Wc.e = ^p,e/10 for the simulation with Vb = 0.6c, and Lo^^e — i^p.e/lO 
and Wc.e = Wp^e/100 for the two simulations with — 0.9c. With these ratios we can 
resolve all relevant electron timescales without requiring long simulation times. The 
weak magnetic field equals that in Ref. 32 and the strong magnetic field that in Ref. 
| 49| . In both these works, the authors considered nonrelativistic beam speeds. 

Typical ratios for Wp^e/wce may be larger in astrophysical environments. We get, 
for example, the ratio Wp.e/wc.e — 3200 if we use the estimated average magnetic field 
B = 10 nT and the estimated average electron number density rie = lO^'^m"'^ for the 
accretion disk orbiting the central black hole in NGC 4261 Spatial and temporal 
fluctuations in the values of B and could, however, lead to lower local characteristic 
frequency ratios in the accretion disk. 

3.1. Slower mildly relativistic beam, vi, = 0.6c 

We resolve the electron distribution by 1600 ppc and each of the proton species by 
784 ppc. The simulation duration is 48GTp and during this time the proton beams 
rotate by an angle a ~ 7.5°. The electric field follows closely that of the unmagnetized 
plasma as shown in figureEl Figure El depicts the time evolution of the kinetic energies 
of the plasma species for unmagnetized and magnetized regimes. The electrons are 
accelerated most rapidly during the initial linear stage of the developing instability 
{t < 50Tp), when they are first trapped by the wave. Here, the electron acceleration is 
almost identical in both magnetised and unmagnetised cases. This is also reflected by 
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Figure 9. Time evolution of the total energy of different particle species in 
units of the initial energy of proton beam 1, Eg, for dj, = 0.6c. The short 
traces correspond to unmagnetized plasma, long traces to magnetized plasma, 
(a) Logarithm of the electron kinetic energies; (b) logarithm of the kinetic energy 
of the background protons; (c) shows the energy of proton beam 1; (d) the energy 
of proton beam 2. 

the similar energy loss of the beams to the background electrons and protons during 
this time. At t « 5QTp the curves diverge. The particle energies of all species in an 
unmagnetized plasma reach a plateau, i.e. the instability has been quenched. Proton 
beam 2, which has been connected to the stronger wave in figure [3 has lost more 
energy. In the magnetized plasma, on the other hand, the electron energy continues 
to rise after t « 50Tp. The electron energy is here a few percent of the beam energy, 
and the continuing acceleration extracts noticeable fractions of the beam energy. The 
ESA extracts energy first from proton beam 2 because the associated strong wave 
traps electrons first. In contrast, the energy of proton beam 1 drops rapidly at 
around t « 250Tp. The energy gain by the background protons is comparable for 
the unmagnetized and magnetized plasmas. 

During the initial stage, beam 2 is dominant in the ESA. The plasma distribution 
at t = 207Tp is shown in figurelTUI The elect rons have reached a spatially isotropic and 
gyrotropic momentum distribution. We discuss below the mechanism by which the 
electrons have been thermalized orthogonally to B, for the simulations with Uf, = 0.9c. 
We find strong momentum oscillations of beam 1 and well-developed phase space holes 
in beam 2. These beam distributions are further evidence of a much more strongly 
nonlinear wave associated with beam 2. 

The rapid drop in the energy of beam 1 is associated with spatial concentration 
leading to the development of a density peak shown in figure llir a) where the local 
beam density exceeds the average density of the background electrons. Strong 
oscillations at 2A„ produce a density modulation which leads to non-linear steepening 
of wave-fronts as they traverse the system. This results in the observed density spike. 
The corresponding "bunching" in phase space indicates the onset of wave breaking 
and is accompanied by a jump in the mean speed of the beam, see figure [TTT c'). The 
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Figure 10. Particle distribution functions at t = 207Tp for Vf, = 0.6c. All 

phase space densities reflected by the colour scale are in units of the logj^Q of 
the number of simulation particles, (a) The spatially homogeneous electron 
momentum distribution; (b) the hot thermalized electron momentum distribution; 
(c) the strongly oscillating proton beam 1; (d) well-developed phase space holes 
in the distribution of proton beam 2. 



development of this density spike is similar to that forming at the leading edge of 
plasma expanding into vacuum. This is detailed in Ref. 34 where it is discussed 
analytically, using a hydrodynamic self-similar and non-relativistic approach, how 
a sharp ion front develops in response to a preceding electron cloud. In such a 
system, the decreasing velocity profile in the expansion direction of a comprcssional 
wave, that expands into a less dense medium implies a steepening and eventually 
breaking of the wave. At the point of wave-breaking, the distribution function becomes 
multi-valued and the fluid approximation breaks down, but a treatment based on the 
Vlasov equation (also detailed in Ref. Pl]) demonstrates the formation of structures 
qualitatively similar to those displayed in Fig lllf c'). Furthermore, beam momentum 
jumps associated with a strong beam charge density modulation, have also been 
reported in Ref. |33|. 

The difference between fjeure iriT c') at t — 250Tp and its earlier counterpart figure 
lior d at t — 207Tp reflects, together with llir al. the key physical process underlying 
localised ESA in the present context. Its consequences for the electron acceleration are 
immediately visible in figure IT^ at t — 25QTp, which is the counterpart of the earlier 
figure EJb) at t — 207Tp. The local density accumulation of beam 1 is connected 
to a strong electrostatic wave potential moving towards increasing values of x. The 
resulting electric field is strong enough to trap electrons, but only over the spatially 
localised high density peak. The result is a localised ESA which leads to an electron 
beam that accelerates out of the heated electron distribution. This electron beam and 
the associated "frying-pan" distribution are shown in figure [T^ 

The local electron acceleration counteracts further proton beam charge 
accumulation and eventually reduces the beam potential, so that the electrons detrap. 
This detrapping arrests both the electron acceleration and the associated energy loss 
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Figure 11. Proton beam distribution functions at t = 250Tp for vi, = 0.6c: (a) 
and (b) show the number density of beam 1 and beam 2, respectively, in the box 
frame of Ref. in units of the background electron number density; (c) and (d) 
show the phase space distributions of beams 1 and 2 respectively. The colour 
scale is in units of the logjg of the number of simulation protons. 




Figure 12. The electron momentum distribution at t = 250Tp for vi, = 0.6c; 
The colour scale is in units of the logj^g of the number of simulation electrons. 
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Figure 13. Proton beam distribution functions at t = 455Tp for vi, = 0.6c: (a) 
and (b) show the number density of beam 1 and beam 2 respectively in the box 
frame of Ref., in units of the electron number density; (c) and (d) show the phase 
space distribution of beams 1 and 2, respectively. Both colour scales are in units 
of the logj^g of the number of computational protons. 



of proton beam 1 in figure|51at t « 270Tp. The beam density peak evolves to the much 
less pronounced but still strongly localized spike shown in figure [TSl At this time we 
also observe the onset of a merger of the two phase space holes in proton beam 2. 

The electrons from the handle of the "frying-pan" distribution shown in figure 
1121 have detrapped at this time, and they gyrate freely in the magnetic field. Their 
relativistic mass increase reduces their rotation frequency for increasing momenta, and 
the accelerated beam electron distribution therefore develops into the spiral shown in 
figure 

3.2. Faster mildly relativistic beam, = 0.9c 

Due to the anticipated larger momentum range accessible to the electrons we use 3200 
ppc for the electrons and 2048 ppc for the protons in our PIC simulations, which 
provides adequate momentum space resolution. We model the plasma for a duration 
of 980Tp. In the simulation with LOp^e/^ce = 10 the proton beams have rotated by 
8.4°. 

The electrostatic field in the weakly magnetic simulation with oip e = lOOt^ce (I'ed 
curves in figure grows similarly to that in the unmagnetized plasma in figure [S] 
The visible rapid fluctuations are similar to those that have been attributed to the 
trapping of electrons [221 ■ The wave in the more strongly magnetised simulation with 
<^p,e = lOti-'ce (blue curves in figure 115(1 does not show these. This is possibly due 
to rapid electron acceleration by the magnetic field. The electrons gain substantial 
energy, and thus relativistic mass, during a bounce period in the wave potential. The 
relativistic mass in this case depends strongly on the initial phase space position of 
the electron, and no coherent oscillations of trapped electrons occur. 

Figure ITHI confirms the more rapid acceleration of electrons in the simulation with 
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Figure 15. Electrostatic waves driven by the proton beams with Vf, = 0.9c. The 
blue trace corresponds to LUp^e = 10iJc,e and the red trace to uip^e = lOOcJce. The 
amplitudes of the waves with negative phase speed {ku < 0) are shown in (a), 
while (b) shows the amplitudes of the waves with ku > 0. 
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Figure 16. The time evolution of the particle energies in units of the initial 
energy of beam 1, Eq, for vi, = 0.9c and for the weak and strong magnetic fields, 
(a) Electron kinetic energies: the upper trace corresponds to the strong magnetic 
field, (b) kinetic energy of the background protons: the strongly oscillating lower 
trace corresponds to the weak magnetic field, (c) Energy of beam 1: the lower 
trace corresponds to the strong magnetic field, (d) Energy of beam 2: the lower 
trace corresponds to the strong magnetic field. 

Up^e = lOwce as follows. The curves representing the kinetic energies of the electrons 
evolve practically identically until t « 125Tp. The initial acceleration mechanism is 
here, as in figure^ the trapping of electrons by the growing wave potential. After this 
time the electrons are accelerated more rapidly in the simulation with ujp^f, ~ lOcJc.e 
and their energy is about one order of magnitude larger than in the simulation with 
Wp^e — lOOcjc.e- As the simulations approach the end time t — 980Tp, both curves 
appear to converge. 

The background protons are only weakly accelerated. The high frequency 
oscillation of the proton energy for cjp e — lOOwce is likely to be the proton reaction 
to the oscillation of the corresponding (red) electric field in figure ^| These high- 
frequency oscillations are almost absent in the case of cjp e — lOwc.e which is also in 
accordance with the absent field fluctuations in figure 1X^1 Both beams lose about 2% 
of their initial energy during their non-linear interaction with the background plasma. 
Of particular interest is the rapid drop of the energy of beam 2 at f « 200tp. 

The energy transfer from the beams for Wp^e = lOOwce accelerates the electrons 
to highly relativistic speeds. In figure [T7I the well-defined phase space structure, with 
distributions that are not yet spatially isotropic, is evidence that electrons are still 
trapped, and therefore undergoing ESA. This explains why electron energy is still 
increasing in figure ITHT a^ at the end of the simulation. 

For [jjp^e = lOwc^e the electron distribution extends to |p| k, ISOrrieC and is spatially 
isotropic at the simulation's end, but not yet gyrotropic, as figure 1181 shows. Both 
beams show phase space holes, which are the reason for the persistent electric field 
amplitudes (blue curves in figure IT^ even at late times. 

The two attached movies show the time evolution of the electron momentum space 
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Figure 17. Phase space distributions for Vf, = 0.9c and ojp^e = 100aJc,e at 
t = 980Tp. (a) Spatially almost isotropic electron distribution reaching a peak 
momentum in excess of |p| = lOOmeC. (b) Electron momentum distribution 
orthogonal to B is shown, (c) Phase space distribution of proton beam 1, (d) 
that of beam 2. 



distribution in the Px,Py plane for Vb = 0.9c and for both magnetic field strengths. 
We observe the electrons along the a;-direction towards increasing values of x\ this has 
the effect of reversing the sense of rotation relative to figures El and ^1 For weak 
magnetic field, Movie 1 shows the ESA for LUp^e = lOOwc.e and Movie 2 shows the ESA 
for Wp,e = 10wc,e. In both movies, the total momentum is mapped into the colour 
scale. Movie 1 indicates that the phase space distribution remains unchanged except 
for a scaling factor that increases approximately linearly as a function of time. Movie 
2 shows that the strong magnetic field case yields a rapid acceleration of electron 
beams which explains the rapid drop of the energy of beam 2 in the lower curve of 
figure ITCT d). It also shows repeated ESA bursts that occur whenever the electric field 
amplitude, shown in blue in figure [T^ rises to a value at which it can trap electrons. 
Another important observation from Movie 2 is the strong increase of the momentum 
spread of the narrow electron beams that are produced by ESA. This beam heating 
is most efficient when the electrons move parallel or antiparallel to the x-axis. In 
this case, they can interact resonantly with the electrostatic waves. This electron 
thermalisation resembles the stochastic acceleration mechanism examined in Ref. j5()| 
for nonrelativistic phase speeds of the waves. 

Finally, we integrate the electron distribution functions in figures [TTT a) and^Ja) 
over space to obtain their phase space density as a function of the total momentum 
IpI, which we show in figure [1^1 

Electrons in the more strongly magnetized plasma reach a higher peak 
momentum. The heating in Movie 2 that, we believe, is due to stochastic heating, has 
apparently led to the development of a power law distribution at high \p\. In contrast 
the momentum distribution reached by the electrons in the more weakly magnetized 
plasma (red trace in figure IT?^ shows an abrupt density drop at high \p\. Here, the 
fastest electrons are still trapped and have a well-defined peak momentum. 
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Figure 18. The phase space distribution for vi, = 0.9c and Ljp^e = 10a;c,e at 
t = 980Tp: In (a) we find a spatially isotropic electron distribution reaching 
a peak momentum \p\ = ISOnieC. In (b) the electron momentum distribution 
orthogonal to B is shown. In (c) we show the phase space distribution of proton 
beam 1, and in (d) that of beam 2. 




Figure 19. The electron momentum density for beam speed vi, = 0.9c at 
t = 980Tp: The blue curve corresponds to the simulation with ujp^e = 10LJc,e, and 
the red to ujp^e = 100aJc,e- The density is given in units of simulation electrons. 
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4. Discussion 

In this paper, we have examined the efficiency of electron surfing acceleration in 
producing relativistic electrons in plasmas. We have used both Vlasov and PIC 
computational approaches. The chosen initial conditions may be representative of 
the foreshock region of plasma shocks in the accretion disks of microquasars, and of 
the central black holes of active galactic nuclei. We have found that across the beam 
velocity band between vi, = 0.6c and vi, = 0.9c, the peak momentum reached by the 
electrons increases by almost an order of magnitude. The most important variable is 
the increasing lifetime of the saturated electrostatic wave. This is also been observed 
in both the PIC and Vlasov codes, which yield quantitatively similar results. This 
similarity, despite the different methods by which both codes solve the relativistic 
Vlasov-Maxwell equations, is evidence for a physical and not numerical origin for the 
increasing lifetime of the saturated electrostatic wave. This lifetime is restricted for 
fast but non-relativistic phase speeds of the electrostatic waves j^|321- In addition, 
the increasing phase speed of the wave, and thus the increasing Lorentz force that 
acts on the electrons, also contributes to stronger acceleration. The relation between 
the peak momentum reached by the electrons and the beam speed is thus a function 
that grows faster than linear. This trend persists up to Vb = 0.99c, for which the peak 
electron momentum has been shown to increase by another order of magnitude jlOj . 
ESA cannot however, accelerate electrons to a peak momentum exceeding mp/rrie |10j . 

We have compared the final momentum distributions for both beam speeds, for 
unmagnetized and magnetized plasmas. The introduction of a magnetic field with the 
ratios Wp^e/'^c,e = 100 and 10 has increased the peak momentum reached by electrons, 
compared to the case of unmagnetized plasma, from |p| /mgC w 2.5 to |p| /mgC w 10 
for Vb = 0.6c and from \p\/meC « 15 to |p| /mgC « 200 for Vb = 0.9c. 

We conclude that ESA works best if, first, the electrostatic waves have a 
relativistic phase speed and, second, the magnetic field is strong enough to accelerate 
the electrons significantly during the lifetime of the saturated wave. In addition our 
simulations show that the overall electron acceleration can be enhanced by secondary 
instabilities that are triggered as a reaction to the strong electron acceleration. 
We have found that the strong wave fields can lead to the local accumulation of 
beam protons which can, by their associated electric fields, yield further localised 
acceleration of electrons. Such an accumulation of the beam density has, in other 
contexts, been investigated analytically for the case of a hydrodynamic plasma, 
expanding into a vacuum, 1841 and for a Q-machine experiment . and it may 
play an important part in electron acceleration to high energies in astrophysical 
environments. This spatially localized acceleration is responsible for the observed 
development of a frying-pan distribution that allows some electrons to double their 
momentum compared to that reached after their interaction with the Buncman wave. 
An inverse power law distribution can arise for the most strongly accelerated electrons. 

The electron acceleration times, which arc of the order of a few hundred plasma 
periods Tp, can yield relativistic electrons from a thermal pool in short times. ESA is 
thus a potential mechanism that could give rise to the rapid fluctuations in the emission 
intensity of radio synchrotron radiation from the accretion disks of microquasars. As 
an illustration we can take the average density of 10^^ /m? reported for electrons in the 
accretion disk of the AGN NGC 4261 which leads to Tp k. lO'^s. The low average 
magnetic field in the accretion disk gives rise to a weak ESA capable of producing 
relativistic electrons from the thermal background over an acceleration time of a few 
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thousand Tp, i.e. a few milliseconds. The rising time of synchrotron emissions would 
then be of the same order. Such a mechanism may explain the rapid fluctuations in 
the synchrotron emission from accretion disks of microquasars jll| or from GRBs l2|. 
However, we have to point out that our plasma initial conditions are idealized and, 
in addition, only limited data is currently available about typical plasma parameters 
of microquasars accretion disks or gamma ray bursts. The present work constitutes 
a first step towards a kinetic model of electron acceleration in the foreshock region of 
mildly relativistic astrophysical shocks. It shows that ESA is a promising candidate 
for the rapid generation of highly relativistic electrons in astrophysical environments 
that support such shocks. Future work should consider ESA driven by ion beams 
that have evolved self-consistently out of a shock, together with the effects of oblique 
magnetic field angles. 
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